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Abstract 

We present our final results for the excited charmonium spectrum from a 
quenched calculation using a fully relativistic anisotropic lattice QCD action. 
A detailed excited charmonium spectrum is obtained, including both the ex- 
otic hybrids (with J^"-^ = 1 '",0^ ,2^ ) and orbitally excited mesons (with 
orbital angular momentum up to 3). Using three different lattice spacings 
(0.197, 0.131, and 0.092 fm), we perform a continuum extrapolation of the 
spectrum. We convert our results in lattice units to physical values using 
lattice scales set by the ^Pi — IS splitting. The lowest lying exotic hybrid 
1"+ lies at 4.428(41) GeV, slightly above the D**D (S+P wave) threshold 
of 4.287 GeV. Another two exotic hybrids 0^ and 2"^ are determined to 
be 4.70(17) GeV and 4.895(88) GeV, respectively. Our finite volume analysis 
confirms that our lattices are large enough to accommodate all the excited 
states reported here. 

I. INTRODUCTION 

The hadron spectrum is one of the most prominent non-perturbative consequences of 
QCD. Over the past two decades, enormous progress has been achieved in developing the 
non-perturbative techniques needed to predict these hadron masses from first principles. It 
is crucial that our theoretical techniques reproduce the observed hadron spectrum to justify 
predictions for more complicated quantities such as the weak matrix elements. 

On the experimental front, data for the excited charmonium spectrum above the DD 
threshold are scarce, and no heavy hybrid signal has been observed yet, even though some 
potential candidates for light exotic hybrids have been found and discussed extensively 
[1]. However, the situation will change dramatically in the next few years with more data 
from the 5-factories PEP-II (Babar) and KEKB (Belle). Even more important results are 
promised by the focused charm physics programs proposed for the upgraded CLEO-c detec- 
tor [2] (expected to start running in early 2003) and the new BESIII detector [3] (expected 
in 2005-2006) at BEPC (Beijing Electron-Position Collider). Both CLEO-c and BESIII 
will collect large amounts of data for the charmonium spectrum above the DD threshold 
and provide badly needed, accurate results for the excited charmonium spectrum, including 
radially and orbitally excited conventional cc mesons and gluon-rich states (glueballs and 
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hybrids). Better theoretical predictions, especiaUy those from lattice QCD calculations, will 
provide important guidance for these experimental efforts. 

In contrast to the light hadron spectrum, the heavy quark spectrum is much better 
understood theoretically with various phenomenological models such as the static potential 
models and, more reliably, with the low energy effective theories such as non-relativistic 
QCD (NRQCD) [4,5]. NRQCD has been the dominant approach to the study of heavy 
quark physics for many years and has been quite successful, especially for the bottomonium 
systems. However, it is rather difficult to control the systematic errors of NRQCD, including 
the relativistic corrections and the radiative corrections, which were shown to be quite large 
for charmonium [6,7], and still sizable for bottomonium [8,9]. Of even greater concern are 
the finite lattice spacing artifacts that could not be controlled by taking the continuum limit 
because of the restriction: M^a > 1. Nonetheless, all these approximations to QCD provide 
important guidance and calibration for fully relativistic lattice QCD calculations and are 
quite useful for the understanding of systematic errors such as those arising from quenching 
and finite volume effects. 

Our ultimate goal is to determine the heavy quark spectrum (extrapolated to the con- 
tinuum limit) from fully relativistic lattice QCD. However, the large separation of energy 
scales {uigV^ -C rUgV -C rUq) in the heavy quark system makes it prohibitively expensive to 
study heavy quarks on a conventional isotropic lattice. Earlier relativistic studies employed 
the "Fermilab approach" [10,11], which uses a space-time asymmetric, 0{a) improved action 
to control the O^Mgo) errors on an isotropic lattice with equal lattice spacings in the space 
and time direction {at = a^). In this study, we employ a newer approach called relativistic 
anisotropic lattice QCD, which is a generalization of the "Fermilab action" to an anisotropic 
lattice with at < a^. The anisotropic lattice gauge action was developed in the 1980s [12], 
and has been used quite successfully for glueball spectrum calculations [13,14] and thermo- 
dynamic studies [15]. However, it was not until the late 1990s that an improved relativistic 
anisotropic fermion action was implemented for Wilson-type quarks [16,17]. Since then, the 
relativistic anisotropic lattice technique has been applied to calculate the heavy quark spec- 
trum with unprecedented accuracy, including both charmonium [17-19] and bottomonium 
states [20]. 

In addition to conventional mesons made of a quark-antiquark pair (within the naive 
quark model), hybrid mesons containing valence gluons are also believed to be present in 
the nature ^. Phenomenological studies using the Born-Oppenheimer approximation [21], the 
bag model [22,23], the flux-tube model [24], and the QCD sum rules [25,26], have provided 
the theoretical evidence and given semi-quantitative estimates for these hybrid states, which 
are further confirmed by NRQCD [27] and standard lattice QCD calculations [28]. In the 
naive quark model, not all J^'-'' quantum numbers are available, since the parity (P) and 
charge conjugation (C) of a conventional qq meson are given by: 

p= (-1)^+1, C= (-1)^+^, (1.1) 

where L and S are the orbital angular momentum and spin of the qq pair. Hybrid states 



We use qq to represent a conventional meson (a pure quarkonium) and qqg for a hybrid meson. 
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possessing so-called exotic quantum numbers (eg. 1 '", O"*"", and 2"*"") have attracted con- 
siderable attention because such states do not mix with the conventional qq mesons and 
therefore should be more easily distinguished in experiments than can non-exotic hybrids. 
More comprehensive reviews of hybrid excitations can be found in Ref. [1,29]. 

However, hke all the other highly excited states, these hybrid mesons are difficult to 
measure accurately on a lattice because their signals are much nosier than those of the low 
lying mesons. Thus, it is crucial to have fine temporal resolution as well as highly optimized 
operators to collect as much information as possible from meson correlation functions even 
at short time seperation. The anisotropic lattice technique provides an efficient framework 
for such high precision calculations, because it allows independent control of the temporal 
and spacial lattice spacings, thereby permitting simulations with a very fine grid in the 
temporal direction. 

In this work, we employ well optimized extended operators to study the excited char- 
monium spectrum, including the orbitally excited mesons (P-, D-, and F-wave) and hybrid 
mesons with exotic quantum numbers. This complements our group's previous charmonium 
calculations on anisotropic lattices [18], which used only local operators and focused on the 
conventional low lying S-wave and P-wave cc mesons. This paper is a follow-up and complete 
analysis of our earlier work [30] . 

This paper is organized as follows. Sec. II describes the anisotropic lattice gauge action 
and fermion action. A discussion of the construction and optimization of meson operators 
is given in section III. Sec. IV contains the details of our simulations, including simulation 
parameters and some measures taken to improve computational efficiency. In section V, 
we present our spectrum results and discuss the systematic errors and the theoretical and 
experimental implications. Section VI contains concluding remarks. In the appendix, we 
describe the notation used in this paper, our meson naming convention, and some relevant 
properties of the hyper-cubic symmetry group. 



We employ an anisotropic gluon action [12], which is accurate up to C(a^, a^) discretiza- 
tion errors: 



This is the standard Wilson action written in terms of simple plaquettes, P^u{x). Here 
labels the time component and the index i = 1, 2, 3 runs over the three spatial directions. The 
parameters j3 and are the bare coupling and bare anisotropy respectively, which determine 
the spatial lattice spacing, a^, and the renormalized anisotropy, ^, of the quenched lattice. 
The continuum limit should be taken at fixed anisotropy, ^ = Cbs/O't- For the heavy quark 
propagation in the gluon background we used the "anisotropic clover" formulation as first 
described in [17]. The discretized form of the continuum Dirac operator, Q = mq + 1^, reads 



II. ANISOTROPIC LATTICE ACTION 




(2.1) 



= - (a^/2)7^A^ . 



(2.2) 
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Here V^^ is the symmetric lattice derivative: 

q{x) a^V^ q{x) = q{x) U^{x) q{x + /i) - q{x) q{x - fi). (2.3) 

and is defined by a'j^A^q{x) = U^{x) q{x + fi) — 2 q{x) + U^^{x) q{x — fi). For the field 
tensor F^^, we choose the traceless cloverleaf definition which sums the 4 plaquettes centered 
at point X in the {/i, v) plane: 

a^a, F^,{x) = I [P^y{x) + P,^(x) + Pp^(x) + P^^(x) - h.c] , (2.4) 

where /2 indicates the negative direction —fi, so Ppp{x) represents W{x — ji,fj,)W{x — fi — 
i), v)U{x — /t — C',^)U{x — z/, u). This is indeed the most general anisotropic quark action 
including all operators to dimension 5 up to 0(3) symmetric field redefinition. We have 
chosen Wilson's combination, W^, of first and second derivative terms so as to both remove 
all doublers and ensure the full projection property. The five parameters in Eq. (2.2) are 
all related to the quark mass, rUg, and the gauge coupling as they appear in the continuum 
action. By tuning them appropriately we can remove all 0{a) errors and re-establish space- 
time axis-interchange symmetry for long-distance physics. A more detailed discussion of 
these parameters can be found in Refs [17,18,20,31]. 

Classical values for these parameters have been given in Ref. [18]: 



mo = mg(l + i asTUq) (2.5) 
1 + i asiriq 

= ^s-—i 2.6 

1 + i atuiq 

Cs = i^s (2.7) 

Ct^Uiys + r^t-) . (2.8) 



2 V ^ a. 

Simple field rescaling enables us to set one of the above coefficients at will. For convenience 
we fix z/g = 1 and adjust Ut non-perturbatively requiring that the mesons obey a relativistic 
dispersion relation ( c(0) = 1 ): 

£;2(p) =£;2(0) + c2(p) p2 + C»(p4)... . (2.9) 

We also choose mo non-perturbatively, such that the spin average of the lowest S-wave 
mesons (IS) matches its experimental value: |M(^S'i) + |M(^S'o) = 3.067 GeV for char- 
monium. For the clover coefficients (c^, q) we take their classical estimates from Eqs. (2.7) 
and (2.8) and augment them by tadpole improvement. 

Cs ^ Cs/{ul) (2.10) 
ct ^ ct/{utul) (2.11) 

The tadpole coefficients have been determined from the average link in Landau gauge: = 
1/3 (tr (a;)) Landau- For brevity we will refer to this as the Landau scheme. Any other 
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choice for (cs, q) will have the same continuum limit, but with this prescription we expect 
only small (^(a^a) discretization errors. 

Our action is a generalization of the "Fermilab action" [10], which can be considered 
as a special case of our action with ^ = 1. It has been demonstrated [17] that the mass 
dependence of the input parameters is much smaller for an anisotropic lattice ( ^ > 1 ), 
which makes the tuning of the parameters easier in the anisotropic case. 

III. MESON MASS MEASUREMENT 

We extract meson masses from the time dependence of Euclidean-space two-point corre- 
lation functions, which is the standard procedure to determine hadron masses from lattice 
calculations. For the hadron source and sink, we use both local and extended meson opera- 
tors constructed from fundamental bilinears in the form of 

M{x) = ^(x) r, V,[U] Vk[U] ^{x) , (3.1) 

here Fj is one of the 16 Dirac F-matrices and Vj is the symmetric spatial lattice derivative 
operator (see Eq. 2.3). To obtain the dispersion relation, -E(p), we project the meson 
correlator onto several different non-zero momenta by inserting the appropriate phase factors, 
exp (— ipx), at the sink, and sum the sink over a spatial hyper-plane: 

M(i,p) =5^M(i,x)exp(-zpx), (3.2) 

X 

where 

p = n with n = (0, 0, 0), (1, 0, 0), (1, 1, 0), (2, 0, 0), (2, 2, 0), . . . (3.3) 

To vary the overlap of a meson operator with the ground state and excited states, we 
also implement a combination of various iterative smearing prescriptions for the quark fields 
and gauge links, which is described below. 

For the quark fields, we use box sources to approximate the spatial distribution of the 
quark field. More specifically, we set the source to be unity within a spatial hyper-cube 
and zero elsewhere. This allows us to vary the extent in a rather straight forward man- 
ner. As such a formulation is not gauge invariant, we have to fix the gauge, and we choose 
the Coloumb gauge for convenience. Alternatively, we can use a gauge invariant formu- 
lation, where repreated action of the gauge invariant Laplacian generates a Gaussian-like 
distribution and provides finer control of the overlap with different states: 

^!{x) ^ *(")(2^) = (1 - -V'[[/])"*(x) (3.4) 

This scheme is referred to as Jacobi smearing. 

For the gauge field, we apply the well established iterative APE-style fink fuzzing algo- 
rithm [32]: 
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staples 

U'^i^) = Tsuis) X^{x) , (3.5) 

where a is the hnk fuzzing parameter, S{x, jj) is one of the six staples around the hnk C/^(x), 
and Vsuia) is a projection back to SU{3) by requiring the new hnk C/^(a;) maximize: 

ReTr {u'^{x)X,{x)) . (3.6) 

The maximization is achieved by several Cabibbo-Marinari pscudo-heatbath update steps 
(fuzzing CM hits). We repeat the fuzzing step defined by Eq. 3.5 for N times (N is the 
fuzzing level). Link fuzzing is quite effective in suppressing gluon excitations and has been 
widely used in studies of glueball and hybrid hadrons. 

This setup allows us to extract reliably both the ground state energies and their exci- 
tations from correlated multi-state fits to several smeared correlators, C*(t), with the same 

nflt 

C%t, p) ^ (M(t, p) M^(0, p)) = J2 «^ (P) (e-^'^''^ * + e-^'(P)(^'-*)) . (3.7) 

i=l 

Since we are working in a relativistic setting, the second term takes into account the back- 
ward propagating piece from the temporal boundary. 

We construct meson operators that have definite lattice quantum number i?^*^, in which 
R is one of the five irreducible representations of the hyper-cubic group: Ai, A2, E, Ti, 
and T2. A brief discussion of the hyper-cubic group is given in appendix C, including 
our choice of bases that carry the irreducible representations and some relevant Clebsch- 
Gordan coefficients. However, determining the continuum quantum number J^'^ of an 
operator is complicated by the fact that the mapping from the finite number of irreducible 
representations of hyper-cubic group to the infinite number of irreducible representations of 
the continuous rotation group is non-unique. Eq. 3.8 summarizes the mapping from R to 
the first few J numbers (up to 4). 



Ai - 


^ J = 


0,4 


A2 - 


J = 


3 


Ti - 


^ J = 


1,3,4 


T2 - 


J = 


2,3,4 


E - 


J = 


2,4 



In most cases, we are only interested in the lowest lying state, which can be extracted from 
a one-state fit at large time separation or from a multi-state fit at shorter time separation. 
However, the lowest state does not always correspond to the lowest J number, the nature of 
a state (orbital excitation, radial excitation, gluonic excitaion, and etc.) must also be taken 
into account. A notable example is meson operator h\ x D_T2 (detailed in Table VII), which 
projects to both J = 2 and J = 3 according to table 3.8. The first guess is that the ground 
state has J — 2 and is degenerate with the ground state of operator 61 x D_E. However, 2+~ 
is an exotic hybrid meson , which is much heavier than the non-exotic 3"*"" state. Therefore, 
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the ground state of bi x D_T2 should have quantum number J = 3 and be degenerate with 
bi X D_A2. which is confirmed by our calculation. Table VII summarizes all the meson 
operators used in this work and the quantum numbers (both lattice and continuum) of the 
ground states. Our naming convention is explained in appendix B2. 

Measuring meson correlators using extended operators is much more expensive than 
those Tising local operators, because more quark propagators arc needed. We take some 
steps to improve the computational efficiency, which are explained in detail in section IV. 

IV. SIMULATION DETAILS 

The strategy and first results for charmonium have already been presented in [17,18]. 
The basic idea is to control large lattice spacing artifacts from the heavy quark mass by 
adjusting the temporal lattice spacing, at, so that m^Ot < 1. We use ^ = 2 lattices with /3 = 
5.7, 5.9, and 6.1. Table I summarizes the parameters of our simulations. 

For the generation of quenched gauge field configurations we employ a standard heat- 
bath algorithm as it is also used for isotropic lattices. The only necessary modification 
for simulating an anisotropic plaquette action is to rescale all temporal links with the bare 
anisotropy, ^o- Depending on the gauge coupling and lattice size, we measure hadron propa- 
gators every 100-400 sweeps in the update process, which is sufficiently long for the lattices 
to decorrelate. 

The renormalized anisotropy is related to the bare parameters through: 

e = (4.1) 

A convenient parametrization for rj is given in Ref. [17] 

/.^ /-. In ^ 1 — 0.55055(7^ , , ^, 

,(e,/3) = l + (l--)„K)^-^^^|. (4.2) 

where fii{C,) is also determined non-perturbatively in Ref. [17] and 77(^,00) = 1 as it ought 
to be. We use this form to determine the appropriate ^0 for each value (3 — Q/g^ — 5.7, 
5.9, and 6.1 for ^ = 2. 

Although the fine temporal resolution of an anisotropic lattice increases the number 
of measurements at short time separations, optimization of the meson operators is also 
quite important. The correlation between measurements of correlation function at different 
time slices reduces the amount of "useful" information for a correlated fit (in the case of 
strong correlation between neighboring time slices, little can be gained by using a higher 
anisotropy) . Using multiple operators (with different optimization parameters) and applying 
multi-state fits can provide stronger constraints for the fits. A source of added difficulty is 
the possibihty that the short distance correlators may be dominated by excited state signals, 
which can generate a false plateau, causing incorrect identification of an excited state as the 
ground state. We observe that the hybrid states are rather insensitive to the box size or 
quark field smearing, but are highly sensitive to gauge link fuzzing. Appropriate link fuzzing 
suppresses excited states because fuzzing reduces short-distance fluctuations of the gauge 
fleld background. This should be particularly important for hybrid states which are expected 
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to have a non-trivially excited gluonic component. We use multiple link fuzzing with the 
same quark smearing for hybrid states. In some cases, we employ fuzzing at the sink only 
to save computational time, tests are done to ensure good coupling to ground states. 

Appendix B 1 describe in detail the calculation of meson correlators. First of all, we can 
reduce the number of quark propagators by reusing several basic quark propagators, since all 
the lattice derivative operators we use are linear combinations of basic derivative operators: 
unit(/), first derivative(V.t(i = x,y,z)), symmetric second derivative {Di = SijkVj'Vk), 
asymmetric second derivative(i?j = eijk'Vj'Vk), and symmetric diagonal second derivative 
can be shown that we can construct all the meson correlation functions 
using these basic quark propagators. Furthermore, we can use a linear combination of 
these basic derivative operators to create a "white noise" source, which creates more than 
one quantum number, and rely on the sink operator to project to the quantum number 
of interest. Of course, using such "white noise" source saves computation time at the 
expense of an increased noise level. We carefully tested various combinations and chose 
those combinations that provide the best computational efficiency. 

V. RESULTS AND DISCUSSIONS 

We present our charmonium spectrum in Figure 1. A clear ordering of meson states 
according to orbital angular momentum and gluonic excitation is shown in Fig. 1, as well 
as the much smaller spin splittings. The low lying S-wave and P-wave mesons lie below the 
DD threshold, while the D-wave and F-wave mesons are above the DD threshold. Even 
higher hybrid excitations (with exotic quantum numbers J^^ = l""*", 0+", 2+~) are found 
above the D**D threshold. A portion of the spectrum (the low lying ^5*0, ^Si, and ■'^Pi) is 
taken from Columbia group's previous work (Ref. [18]). 

The lattice scale is set by the ^Pi — 15* splitting ^. Alternatively, we can use the Sommer 
scale ro [33], which was employed in Ref. [18]. Lattice scales from both methods (listed 
in Table I) agree well (within 2-4%) for charmonium, which, perhaps is not too suprising 
since tq is fixed from a potential model whose parameters were chosen specifically to repro- 
duce the experimental charmonium spectrum. However, the discrepancy is rather large for 
bottomonium. We choose the ^Pi — IS scale for this study to be consistent with previous 
NRQCD calculations. 

A. Extrapolation to the continuum limit 

The continuum limit has to be taken to remove the lattice discretization errors. This is 
also a major advantage of relativistic over non-relativistic QCD approach. 



^We use ^Pi instead of the spin averaged P-wave meson mass IP = l/9('^Po + 3^Pi + 5^P2) 
because ^P2 requires non-local operator which is much nosier and ^Pi is extremely close to IP ( 
iPi - IP = 0.86 MeV < IP - IS) 
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The continuum extrapolation of higher excitations is carried out as follows. We first 
calculate the dimensionless ratio of splittings between an excited state and the IS* state (the 
spin-averaged mass of S-wave mesons): R{as) = (M — M{lS))/{M(^Pi) — M{1S)) for each 
lattice spacing a^, and then extrapolate the ratio R{as) to the continuum limit = 0. In 
the end, the continuum result of the excited state is obtained by: 

M\a,=0 = RL=0*CPl-lS)exp + M{lS)exp (5.1) 

where (^Pi - lS)exp = 458.2 MeV and M(lS)exp = 3067 MeV are experimental values. 
Because the splitting between an excited state and the IS" state is rather insensitive to the 15" 
mass (or the quark mass), the imperfect tuning of IS* to its experimental value has negligible 
effect on the physical values obtained from Equation 5.1. However, the incorrect running 
of the coupling in a quenched simulation will still cause a mismatch of scales perceived by 
excited states and by the ^Pi — IS splitting. In contrast, the fine structure sphttings are 
very sensitive to the quark mass, so the deviation of the IS* mass from its experimental value 
should be taken into consideration when quoting the physical values for the spin splittings. 

Using tadpole improved coefficients, we are able to remove the leading 0{as) errors, which 
has be shown by previous studies [17,18,34]. Therefore, we extrapolate to the continuum 
limit assuming that the leading lattice spacing error is 0{al). 

B. Spin splittings 

We first discuss our results for the spin splittings to demonstrate the necessity of a fully 
relativistic treatment, which is crucial for studying the fine structure of charmonium. 

The Columbia group has studied in detail this spin sphttings [18], but the state 
was missing because only local operators were used. Here, we measure all three P-wave 
triplet states ^Pj (J = 0, 1, 2) and calculate the P-wave triplet splitting ratio Rfg = (^P2 — ^ 
Pi)/(^Pi — ^ Pq), which is highly sensitive to relativistic corrections but less affected by the 
systematic errors originating from fixing the overall scale in a quenched calculation. It has 
been shown that the velocity {v /c) expansion of NRQCD converges poorly for charmonium 
spin splittings [6,7]. Our Rfg results are given in table HI. The result at /3 = 6.1 has a large 
error because the plateau of the fit sets in at a rather large temporal separation. Our result 
for the P-wave triplet splitting ratio Rfs is 0.47(13) in the continuum limit, which agrees 
well with the experimental value of 0.478(5). A comparison with a NRQCD calculation [35] 
is given in Figure 2. Our result is also consistent with CP-PACS result from anisotropic 
lattices with anisotropy ^ = 3 [34]. For a detailed discussion on the hyperfine sphtting, we 
refer the reader to Ref. [18]. 

C. Excited states 

We first discuss our result for exotic hybrid states, in particular the 1^+ meson. Figure 3 
is an effective mass plot of the exotic hybrid meson l""*", with each data set corresponding to a 
particular fuzzing parameter a as defined in Eq. 3.5. With a very fine temporal resolution. 
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wc arc able to obtain many measurements at short time separations to facilitate multi- 
state fits and extract both the ground state and excited state masses more precisely. By 
employing different quark smearing and gauge link fuzzing parameters, we obtain different 
couphng to the ground and excited states (referred to as different "channels" later). With 
multiple channels, the multi-state fits are quite stable. We fit the data to 2-state and in 
some cases 3-state ansatz to obtain the ground state mass. Figure 4 shows the multi-state 
fits for meson 1^+. the ground state mass obtained from 1-state fit at large time separation 
{t > 7) is consistent with the results from a 2-state fit {t > 4) and a 3-state fit {t > 1), 
which rules out contamination of the ground state by excited states and clearly demonstrates 
the advantage of anisotropic lattices. For the hybrid states, we did not observe as strong 
a dependence on the quark field smearing as seen for the conventional qq mesons, which 
indicates that the excitation is likely to originate from the excitation of gluonic constituents. 
For a self-consistency check, we measured the 1 ^ exotic meson with 3 different operators 
(oq X V, 6i X V_Ti, and p x B_Ti) and obtained consistent values. Fig. 5 shows the 
continuum extrapolation {0{a^)) using three different operators, with the continuum value 
being 4.428(41) GeV, 4.409(75) GeV and 4.41(18) GeV for p x B_Ti, a'^ x V, and bi x V-Ti 
irrespectively. The operator p x B_Ti has the best signal and smallest statistical error, so 
we quote the final result 4.428(41) GeV for the l""*" mass from this operator. 

Various models of QCD, such as the bag models [22,23], the flux-tube model [24], and 
the QCD sum rules [25,26] have been used to study these hybrid states. For the charmo- 
nium hybrid meson, the adiabatic bag model predicts a 1^"^ ccg state of 3.9 GeV [23] and 
the QCD sum rules predicts 4.1 GeV [26]. We show the comparison of our 1""*" result to 
previous lattice results [1] in Fig. 7. The result from isotropic lattices (MILC [28]) has 
large systematic uncertainties due to a possible contamination from excited states. Previous 
results from NRQCD calculations [27,36] are quite close to ours, which can be explained 
by the expectation that the quarks in a hybrid meson move more slowly than those in the 
lowest lying cc states (as indicated by the hybrid potential study [37]), thus making the 
relativistic corrections small. However, the NRQCD calculations are done on rather coarse 
lattices and the masses can not be extrapolated to the continuum limit to control the lattice 
discretization errors. 

In addition to the exotic 1"+ meson, we also measured two other exotic hybrid mesons: 
O"*"" and 2~^^ . They are heavier than the 1^+ meson, with masses of 4.70(17) GeV and 
4.895(88) GeV, respectively. No clear signal is obtained for , which is predicted by 
phenomenological models to be even more massive than the 2+^ state. The QCD Sum 
rules gives an estimate of 5.9 GeV [38] for the state, around 1 Gev above the 2+~ state. 
Clearly it is challenging to extract such a huge excitation from the noise. Higher anisotropics 
may help, but a more efficient operator than our ai x V-Ai is probably also necessary. 

Most of our meson operators which can potentially couple to non-exotic hybrid states 
actually mix strongly with the conventional qq mesons with the same quantum number 
J^'^, as indicated by the degenerate masses from different operators shown in Table II. For 
example, operator tt x B_Ti has an explicit chromo-magnetic field, making it a potential 
operator for the hybrid meson 1 . However, our calculation shows that it has the same 
mass as the conventional 1 (see Table II). The only exception we find is the 2^+ signal 
from p X B_T2, which is a very good candidate for a non-exotic hybrid meson. First, its 
dependence on quark field smearing and gauge link fuzzing is identical to the exotic hybrid 
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mesons l""*" from p x B_Ti. Second, its mass is much larger than the conventional qq meson. 
However, further mixing test must be done to rule out the possibility that it is the radial 
excitation of L = 2 conventional 2~+. 

A selection rule analysis [39,40] shows that the width of a charmonium hybrid meson is 
narrow if it lies below the D**D threshold (4.287 GeV). Our result for the exotic hybrid 
is 4.365(47) GeV, slightly above the D**D threshold, so the hybrid meson is expected to be 
broad. However, the conclusion is not final due to the ambiguity of scale setting procedure 
for quenched calculations. 

As a new result from lattice QCD, we also obtain reliable predictions for higher, orbitally 
excited mesons with L — 2 and 3 (D-wave and F-wave), which provide more information 
about the long range interaction than S-wave or P-wave mesons. The D-wavc mesons 
2 ^, 2 , 3 ) are around 3.9 GeV, above the DD threshold (3.7 GeV), and F-wave mesons 
(JPC = 3++^ 3+-) at 4.15 GeV, above the DD threshold, but still below the D**D 
threshold. 

D. Systematic errors 

The excited states have larger spacial extent than the lowest lying states, so the results 
for excited states are more likely to suffer from finite volume effects than the lowest lying 
states. Therefore, we compare the masses of all excited states obtained from two different 
lattice sizes {Ng = 8 and Ng = 16 at P = 5.7, { = 2) and find almost no finite size effect, as 
shown in table IV. This is consistent with previous NRQCD study [36,41], which found no 
finite volume effect for a spacial extent larger than 1.2 fm. So we conclude that there is no 
discernible finite vohimc effect on our excited state masses. 

The second source of error is the imperfect tuning of the bare velocity of light Ut, which 
is found to have a rather big effect on the fine structure [18]. Using a. (3 = 5.7, ^ = 4, 8^ x 64 
lattice, we vary i/f by about 10%, which corresponds to about 5% change in the renormalized 
velocity of light. The resultes are given in table V, the changes in excited masses are all less 
than 2% and negligible within errors. 

The quenched approximation is the largest source of systematic error. The running of 
the coupling constant is not correct without the dynamic light quarks. The ratio R2S = 
{2S — lS)/(}Pi — IS) from our quenched lattices is 1.55(9), which is 15% higher than the 
experimental value of 1.34. It affects our result in two ways. First, we convert our lattice 
results to physical units with lattice scales set by the ^Pi — IS splitting, which is ambiguous 
due to quenched approximation. Second, we tune our quark mass to the charmed quark 
mass by matching the 15* mass to its experimental value of 3.067 GeV using this ^Pi — IS 
scale, which in turn may affect our final results. For our excited spectrum, the later can 
be neglected because the mass splitting between a excited state and lowest S-wave meson 
mass IS* is insensitive to quark mass. To demonstrate this mass insensitiveness, we measure 
the excited states on a /3 = 5.7, = 4 lattice with two different quark masses m = 0.32 
and m = 0.25 (corresponding to IS mass of 3.383 GeV and 3.029 Gev) and find no mass 
dependence (as shown in table V. However, for the spin splitting which is approximately 
inversely proportional to the quark mass, the effect is large (eg. a 5% correction (increase) 
in scale will cause a 5% increase in the l^S" mass, which makes the spin splitting in lattice 
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units 5% smaller than what is should be, so the final result for the spin sphtting is 10% 
larger) . 

VI. CONCLUSION 

We have accurately determined the masses of many new and yet unobserved charmonium 
states (both conventional cc mesons and hybrid ccg mesons) from first principles using a 
fully relativistic anisotropic lattice QCD action. While all the new states are predicted 
above the DD threshold, we believe that the exotic hybrids may still be sufficiently long- 
lived to warrant an experimental search in the energy region near and above the D**D 
threshold of 4.287 GeV. This will be an important effort to validate non-perturbative QCD 
at low energies. We add further evidence that a fully relativistic treatment of the heavy 
quark system is well suited to control the large systematic errors of NRQCD. In addition, 
the anisotropic lattice formulation is a very efficient framework for calculations requiring 
high temporal resolutions, which is crucial for the study of high excitations. The quenched 
approximation, which is the major remaining systematic error in our calculation, will be 
addressed by future full QCD calculation. Full QCD anisotropic lattice action has been 
implemented for staggered fermion and used for thermodynamics study [42]. In addition to 
the heavy-heavy system, an anisotropic lattice is also a natural framework to study heavy- 
light systems. 
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APPENDIX A: CONVENTIONS 

In this section, we describe the conventions used in this paper. The variable x specifies 
the coordinates in the 4-dimensional (space-time) volume, with extent Ng along the spacial 
directions and Nt along the time direction. The spacial physical lattice size are given by 
Lg = * tts- X is the 3-dimensional spacial coordinate, and t is the time coordinate, ^'(x) 

or q{x) represents a fermion field. 

The definitions of Sj^fe and Sajk (ct=l,2) are: 

^ijk — \^ijk\ (-^l) 
Sajk = 0(j ^ k), Sni = 1, ^122 = -1, -5222 = 1, ^233 = "1 (A2) 

APPENDIX B: MESON CORRELATOR AND MESON OPERATORS 

1. Meson correlator calculation 

We describe how to construct meson correlators from quark propagators. For simplic- 
ity, we discuss the correlation function between sink meson operator A4{t) — ^2{^2,t) Oi 
ra^i(xi,i), and source operator M(io) = ^2(y2,^o) Oj r^^i(yi,io), in which and F^ 

arc a combination of Dirac 7-matrices, and Oi and Oj are lattice derivative operators. More 
generally, the source operator can be any operator (or a combination of different meson 
operators) that creates the state of interest. 
The correlator is calculated as follows: 

{M^{t)M{to)) = (*l(xi,i)Ft ol 74^2(x2,i) %{y2:to) Oj F^*i(yi,io)) 

= (*l(Xi, i)74rb4 bl *2(X2, t)^2(y2, ^o) Oj F^*i(yi, ^q)) 

= (74r^74 G2(x2,t;y2,to) ^/s (^(yi, to; xi, t) Oi) 

= (7574rl74 G2(x2,t;y2,io) r^75 [Oi ^(xi , t; yi , to)]'^) 

= (F', G2(x2,i;y2,to) [& G(xi,t;yi,to)]^) (Bl) 

Here, we first set up the source (box or smeared) and apply the derivative operator Oj, 
then we calculate the quark propagator ^(yi, to; Xi, t) — (^i(xi, t) Oj ^i(yi, ^o))- Utilizing 
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G'(yi, to; Xi, = 75G^(xi, yi, io)75, we arrive at equation Bl. Then we apply the sink 
derivative operator Oi, insert gamma matrices T'^ = 7574r]^74 and = Tj^j^, and do 
color and Dirac indices contraction. Wc need two quark propagators: G(xi, t; yi, to) and 
G^2(x2, y2, ^o)- In practice, we can reuse the same quark propagators for multiple meson 
operators. 



2. Meson operators and naming convention 

Meson operators with no derivative operators are named ao, , tt, p, ai, and &i, as 
listed in the second row of tabic VI. A meson operator with a lattice derivative operator is 
named X x Y_R, where X represents the gamma matrix, Y is the derivative operator (V, D 
or B), and R indicates the Oh irreducible representation of the meson operator (R is usually 
omitted if it is unique). 

APPENDIX C: HYPER-CUBIC GROUP 

The five irreducible representations of hyper-cubic group Oh and our choice of basis 
vectors for these representations are: 



A2 
E 

T2 



^,1= (:c' + y' + z') 

■0^ = xyz 

ii^Ul)^ {x'-y\y'-z') 

(V'x>'0^,'0z) = {x,y,z) 



(CI) 



(V'i, -^L' i^L) = (y^ + zy^ + ^y + y^) 



The decomposition of often used direct products of hyper-cubic group irreducible repre- 
sentations: 



(C2) 



Clebsch- Gordon coefficients are defined as: 



V'(J,M)= ^ C{J,M\Jl,ml;J2,m2)u{Jl,ml)v{J2,m2) 



(C3) 



mi,m2 



Table VIII, IX, are Clebsch- Gordon coefficients used to construct the meson operators 
in table VII 
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TABLES 





(5.7, 2) 


(5.9, 2) 


(6.1, 2) 


(5.7, 2) 


/at a r \ 


(8, 32) 


(16, 64) 


(16, 64) 


(16, 32) 


coniigs 


1950 


1080 


1010 


1200 


separation 


100 


200 


400 


100 


at (ro) [GcV] 


1.905(5) 


2.913(9) 


4.110(22) 


1.905(5) 


(^Pi - 15) [GeV] 


1.945(26) 


3.021(34) 


4.292(49) 


1.924(25) 


r r 1 

as [tm\ 


0.1974(26) 


0.1307(15) 


0.0920(11) 


0.2052(27) 


Ls [fm] 


1.624 


2.091 


1.472 


2.283 


Co 


1.654729 


1.690713 


1.718306 


1.654729 




0.7762 


0.8091 


0.8280 


0.7762 




0.9394 


0.9504 


0.9569 


0.9394 


/ 

Uot/Uos 


1.2103 


1.1746 


1.1557 


1.2103 




1.208657 


1.1829329 


1.163937 


1.208657 




U.51 


0.195 


0.05 


0.51 


(z^s,l't) 


(1, 1.01) 


(1, 1.09) 


(1, 1.12) 


(1, 1.01) 


Cs 


2.138 


1.889 


1.7614 


2.138 


Ct 


1.3252 


1.2055 


1.1431 


1.3252 


c(0) 


1.000(2) 


0.984(1) 


0.984(3) 


0.991(3) 


Box size 


7 


16 


16 


7 


Fuzzing a 


(1/30, 1/6, 2/3) 


(1/30, 1/6, 1) 


(1/30, 1/6, 1) 


(1/30, 1/6, 2/3) 


Fuzzing Level 


5 


6 


7 


5 


Fuzzing CM hits 


8 


6 


6 


8 



TABLE L Charmonium spectrum simulation parameters. 
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l'5 9 2") 


(P, 1 2') 


r'r\T\i'iT\T\ 1 1 Til ( ( n'^ W 


^'^\T\i'^T\^^^'\ \ii ( 1/ ( n 1 1 

L/UilUililX Lllil I I Lt 1 1 




jPC 




Fr. (CinV) 


Ro (CipV) 

-C/Q 1 v;rt^ V I 








U 






4 S2l'1 7^1 


4 70^ 7^1 


4 84f97'i 









4 40fi^'fi4'l 




4 40Q('75'l 




t'l A V — L \ 




^.O^l -LU 1 






4 41 CIS^ 

^.^X 1 XO 1 


4 51 ^'34'l 


U y\ 1 J L ]^ 


J. 




4 'i94('97'l 




^.4:Z(0\ t:X 1 




n' X n 

Lin /N J--' 




4 4Q4l'Q7'l 


4 71 7l'54'l 


4 9^1 


4 89l'1 0^1 

'X.OZ; I XU J 


5 oqn Q) 


ni X R To 




4 71 Dl'fiD'l 

4:. 1 XUl VJVJ J 


4 S6Sl'54'l 

t:. QUO 1 Ut: J 


5 03^ 1 

(J .UO 1 X X J 


4 895('88') 


5 06f1 P,) 






4 ^\ (^ A\ 


4 fi74l'7'^'l 


"Xr.OOl XZ( I 


'X. i Ul XZr I 


A Q^^^9'^^ 


TT 


u 


O.UO 1 Ul ( 


^ 1945n0'l 


^ 194Q('91'l 

(J. XZi^tJl X J 


O.UUOZil Xc/ J 


o.uuc/t/i 00 1 


ao 


n++ 




'iQ44(''^1 "l 


3 i=>73n('49'l 

0.01 OUl J 


AiOy-±.0 j 


O.fiOU 1 1 t/X ( 




n++ 

U 






3 'i30('1 9") 


O.'XZiUl X X 1 


0.'4:0Ul XO 1 




1 




Q on'^'^i'l 9"! 






o.uooul 00 ) 


^ ^ R 


1 


Q 1 fi4ni'27'l 


Q 904'^0('fi0'\ 




.\jO\j'-k\AO J 


o.uou ( I 00 ) 


P 


1 — 


3 1fi4fi('2n'l 


3 9053('10'l 


3 1Q80('91 

Kj. ±cf(jyj\A ± J 






n Y D 7"i 


1 


'1 ] -. 7/ r - > \ 


• i 0(1 ir;/ 


907f \X\ 

•).A,) I J 


'{ 1 0i']{''\7) 

'.) . _l_ _■ ..7 I * ) 


'{ '?(']R{>( 1 1 

. A\JO I ± I 


rti V R T, 

CfcX A IJ L \ 


X 


0.0 1 tj I Z X ) 


Q fiztl ('9(^^ 

O.U^X I ) 




Q ^'=iR('9Q'l 

.OOOyAiJ J 


0.000 1 ^0 J 




X 


o.uuuui X I 




o.uuoi xo y 


o.ooooi / J 


o.o^y 1 XO J 


U\ A 2^ L \ 


X 




o.u^uui yx I 


O.UXUl X / J 


0.(JUO\ XO J 


O.OKJOyA ± J 


IT V \/ 

/I A V 


X 


0.0\JA\±\J I 


O.UUtJl XZ( I 


O.UUZit X 1 J 


Q zLQnn 7^ 


0.'4:c/0l ox J 




1++ 

X 


O.Ot/ 1 tJl ux / 


O.W^JWt:! .^X / 


^ fi9fi9('fi'^'l 


^ '^^94(''i9'l 

KJ.tJtjAj'^ltJAj J 


^ "i^i (^ ^\ 

0.00 ± I X X J 


til ^ 1 


1++ 

X 


O.OOiJ ± 1 1 ( 


O.UUZil XO ( 


O.UXUXl OZ( J 




O.^c/ZpI XU J 


fj /\ V — L \ 


1++ 

X 


0.00\JA\A i J 


o.(jyoui yx I 


0.(JOU\ X X J 






n-t X \7 F! 


9 — 

Ai 




Q Q1 zL^l 

0. y X / I J-t: J 


Q qqqM 7^ 


0.0\JO[A\J J 


o.oxu I 1 1 


til A V — ^ 2 


o 

A 


o.oy X I ox I 


0. iyXUl X(J I 


O.OiJO \10j 


Q 779/'! R\ 


0. 1 001 OU J 


n V n To 


9 — 

Ai 


O.Oc/XUl c/O I 


0. c/'XOl XO 1 


o.yziOi xo y 


O.OXOl XO J 


O.OXfil J 


TT V n To 


Ai 


O.OO^Ai\ijO J 


O.e/OX I X^ J 


O.yZfOl XO J 


O.OX^l X^ J 


O.OX / I ZiO J 


/n V R To 


2"+ 

A 


^.O'-tOyA L J 




^.\jAiO\\J^ J 


^.0 i Ul OU J 




(_A 1 A 


2++ 

A 


.\J\J\J A \OiJ J 


o.uoux I 1 ) 


Q fi497/'7*^'\ 


Q '^9'^7i'7S^ 


O.OZOt XZ 1 


ai X D_T2 


2++ 


3.6033(84) 


3.6324(76) 


3.6367(60) 


3.5202(74) 


3.520(13) 


p X v_r2 


2++ 


3.5570(62) 


3.6207(91) 


3.605(12) 


3.508(12) 


3.521(19) 


p X i:)_y42 


3~ 


3.912(12) 


3.9760(91) 


3.958(17) 


3.865(15) 


3.879(27) 


6i X D_A2 


3+- 


1.153(29) 


4.333(39) 


1.271(17) 


4.239(48) 


4.312(79) 


hi X L»_T2 


3+- 


4.157(27) 


4.339(42) 


4.383(39) 


4.326(43) 


4.450(70) 


ai X D_A2 


3++ 


4.151(56) 


4.248(54) 


4.344(99) 


4.215(84) 


4.31(15) 



TABLE II. Charomonium spectrum results expressed in GeV. The meson operators are sorted 
by the continuum quantum number {J^^) of the ground state. Only the ground state masses are 
listed in this table. The lattice scale is set by the ^Pi — IS splitting. Continuum results from 
0(0^) and 0{as) extrapolations are listed in the last two columns. The continuum extrapolation 
procedure is described in section VA. The errors quoted here do not include the uncertainties in 
the physical scale. 
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(AO (5.7,2) (5.9,2) 


(6.1, 2) 


cent. 0{a'^) cent. 0{a) 


Rfs 0.474(50) 0.46(8) 


0.52(20) 


0.47(13) 0.47(21) 


R2S 1.50(12) 1.563(45) 


1.517(77) 


1.55(9) 1.54(18) 


TABLE III. P-wave spin splitting ratio Rfg = (^^2 


— ^ Pi)/(^Pi — ^ Pq) and excited state mass 


ratios (i?^ = {X - lS)/{^Pi - IS)). 






{(3,0 (5.7,2) 




(5.7 ,2) 


Ns 8 




lb 


Ls [fm] 1.624 




3.283 


2-+ 0.739((9) 




0.758(13) 


3++ 1.006(56) 




0.992(31) 


1"+ 1.294(16) 




1.281(15) 


0+- 1.351(51) 




1.337(37) 


2+- 1.564(63) 




1.575(71) 



TABLE IV. Finite volume effects of excited charmonium spectrum. The values quoted are 
masses relative to the \S state. 



(A 6 


(5.7, 4) 


(5.7 ,4) 


(5.7, 4) 




0.32 


0.25 


0.32 




1.37 


1.37 


1.25 


2-+ 


0.801(14) 


0.807(15) 


0.801(12) 


3++ 


1.062(23) 


1.078(43) 


1.071(46) 


1-+ 


1.316(31) 


1.312(31) 


1.318(68) 


0+- 


1.461(46) 


1.475(72) 


1.496(57) 


2+- 


1.692(46) 


1.674(61) 


1.715(63) 



TABLE V. Mass and velocity of light dependence of excited charmonium spectrum. The values 
quoted are masses relative to the \S state. 
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X r 


1 


74 


75 


(7i,747i) 


75 7i 


(74757i>7i7j) 


1 


ao = 0++ 


«o = 0+- 


TT = 0- + 


p = l— 


ai = 1++ 


bi = 1+' 


Vi 


1~ 


1-+ 


1 + - 


(0,1,2)++ 


(0,1,2) — 


(0,1,2)-+ 


Bi 


1+- 


1++ 


1 — 


(0,1,2)-+ 


(0,1,2)+- 


(0,1,2)++ 


Di 


2++ 


2+- 


2-+ 


(1,2,3)— 


(1,2,3)++ 


(1,2,3)+- 



TABLE VI. A quick reference of qq and qqg states (exotics in bold font) accessible with 
O X r construction of meson operator, O is a derivative operator and T is a combination of Dirac 
7-matrices. Here is the first derivative operator, = SijkVjVk (i=l,2,3), and Bi = eijkVjVk, 
where Sijk is the symmetric tensor defined by Eq. Al. 
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u 


riA V V7 4 1 






-'2 


— 


CXI V — £ 2 
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o 
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(X\ A U 
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/^/i^ £1 ■ ■ 1 ■111 

Ibf-ijkljJ-'k 


-'2 


9++ 
z 
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o 
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lUbSajkljDk 
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-t 1 
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1 -^/r- fl ■ ■ 1 1^/ •111 
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Q+- 
O 


u\ A t 2 
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TABLE VII. Meson operators, names and quantum numbers. The quantities Sijk and Sajk are 
defined by equations AI and A2. 



20 





4 4 


4 4 


4 4 


4 4 


4 4 


4 4 


4 4 


4 4 


4 4 




"'x ^x 


"'X 




7/ 7; 




11 1) 


7/ 71 


7/ 7T 

"z fj^ 


7/ 11 


V 


i/ V-J 


n 
u 


n 

U 


n 
u 




n 
u 


n 
u 


n 
u 


i/ V<J 


Wl 


1 / /o 


U 


U 


U 


1 /o 

-1/2 


U 


U 


n 
U 


1 /o 

-1/2 


W2 


1 /o 

1/2 


u 


U 


U 


1 /o 

1/2 


U 


U 


U 


-1/V2 





































-1/^/2 











1/^/2 













1/V2 





-1/V2 

















Tyz 
















1/V2 





1/V2 













1/V2 



















rxy 





1/V2 





1/V2 


















TABLE VIII. Clebsch-Gordon coefficients of Ti Ti = Ti 
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TABLE IX. Clebsch-Gordon coefficients of Ti (g) = ^2 £^ 
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FIG. 1. Quenched charmonium spectrum. The experimental values are shown as short hori- 
zontal lines. The long horizontal dashed lines mark the DD and D**D thresholds. The mesons 
(please refer to appendix B2 for our naming convention) plotted are: 7r(0 p(l ), ), 
p X V_^i(0++), p X V_Ti(l++), p X V_r2(2++), vr x D_T2{2-+), p x D_T2{2—), p x D_^2(3— ), 
ai X D_A2{3++), bi x L>_yl2(3+"), ai x S_Ai(0+-), p x B_Ti{l-+), and oi x S_r2(2+-). Some 
of the low-lying meson (7r(0 p{l ), )) are taken from Coumbia group's previous work 

[18]. The lattice scale is set by the ^Pi — 15 splitting. The numerical values are listed in table II. 
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FIG. 2. P-wave triplet splitting ratio Rjs = {^P2 —'^Pi)/{^Pi — '^^'o)- Comparison is shown 
with NRQCD results with relativistic corrections up to 0{mv^) and 0{mv^) [35]. 
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FIG. 3. Effective mass of the 1 exotic meson for /? = 6.1 and ^ = 2 with different APE-style 
fuzzing coeffients a as defined in equation 3.5. 
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FIG. 4. Masses from multi-exponential fits of p x BJ'i propagator (the ground state is the 
exotic 1 ^) for /3 = 6.1 and = 2. The upper band of solid circles and two added upper bands of 
solid triangles show the effective excited state masses. 
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FIG. 5. Continuum extrapolation (as outlined in section VA) of the exotic hybrid 1 ^ (mass 
above IS*) with different meson operators ( p x B_Ti, q'q x V, and bi x V_Ti). 
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FIG. 6. Continuum extrapolation of the exotic hybrid 0^ (oi x B_Ai) and 2+ (ai x B_T2). 
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FIG. 7. Comparison of charmonium hybrid meson 1 ^ mass (relative to the spin-averaged 
S-wave meson mass IS*) as determined by this work, NRQCD (CP-PACS [36]) , and isotropic 
relativistic lattice QCD (MILC collaboration [28]). 
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